Curso II – SISSA

Introducción a Machine Learning - Parte I

Objetivo y contenido

Objetivo y contenido

Objetivo: trabajar con algoritmos de machine learning en R para la predicción espacio-temporal de variables.

  • Introducción a los algoritmos de machine learning
  • Paquetes de machine learning en R
  • Random Forest
  • Support Vector Machines
  • Redes Neuronales
  • Uso de machine learning para la predicción de variables
  • Ejemplo: machine learning para la predicción del caudal de los ríos

Introducción a los algoritmos de machine learning

Machine Learning

Machine learning es una rama de la Inteligencia Artificial que tiene como objetivo utilizar datos y algoritmos para imitar el aprendizaje de los seres humanos (IBM), y se enfoca en realizartareas conocidas.

Machine Learning

Los algoritmos de machine learning se han utilizado ampliamente en:

  • Asistentes virtuales
  • Vigilancia de video
  • Redes sociales
  • Reconocimiento de patrones
  • Análisis de datos espaciales
  • Evaluación de riesgos

Machine Learning

El machine learning ha demostrado ser útil en muchos campos. ¿Cuáles son los pasos para aplicar métodos de machine learning?

  • Recopilación de datos
  • Preparación y preprocesamiento de datos
  • Selección del modelo o modelos
  • Entrenamiento del modelo o modelos
  • Evaluación del rendimiento
  • Mejora del rendimiento del modelo o modelos

¿Qué algoritmo de machine learning elegir?

No existe un algoritmo de machine learning que sea el mejor en todos los casos. Su rendimiento depende del tipo de datos y los objetivos finales. Básicamente, los dos principales grupos de algoritmos de machine learning son:

Algoritmos de machine learning

Los siguientes algoritmos se utilizan a menudo para el aprendizaje supervisado:

  • Nearest neighbour
  • Random Forest (RF)
  • Support Vector Machine (SVM)
  • Redes Neuronales (NN)
  • Naive Bayes

Algoritmos de machine learning

Los siguientes algoritmos se utilizan a menudo para el aprendizaje no supervisado:

  • k-nearest neighbour (K-NN)
  • Análisis de Componentes Principales
  • Descomposición de Valores Singulares
  • Análisis de Componentes Independientes

Paquetes de machine learning en R

Paquetes de machine learning en R

En este curso, nos centraremos principalmente en algoritmos de aprendizaje supervisado, ya que queremos trabajar con regresión y clasificación de variables.

  • RandomForest: implementación de RF
  • e1071: clustering difuso, SVM, Naive Bayes y K-NN
  • kernlab: métodos de machine learning, SVM, clustering espectral y PCA
  • caret: conjunto de funciones que buscan agilizar el proceso de creación de modelos predictivos

Paquetes de machine learning en R

  • DataExplorer: proceso automatizado de exploración de datos, para que los usuarios puedan centrarse en comprender los datos y extraer ideas
  • superml: proporciona una interfaz estándar para usuarios que usan R y Python para construir modelos de machine learning
  • mlr3: programación eficiente y orientada a objetos en los componentes fundamentales del machine learning
  • party: árboles de inferencia condicional

Paquetes de machine learning en R

  • rpart: particionamiento recursivo para árboles de clasificación, regresión y supervivencia
  • neuralnet: redes neuronales utilizando retropropagación, retropropagación resiliente
  • ranger: una implementación rápida de Random Forest, especialmente adecuada para datos de alta dimensionalidad

Random Forest

Random Forest

Random Forest (RF; Biau and Scornet, 2016, Breiman, 2001, Prasad et al., 2006) es un método que se puede utilizar para tareas de clasificación y regresión supervisadas construyendo árboles de decisión utilizando la relación entre variables independientes y dependientes.

  • Es preciso y capaz de manejar muestras pequeñas y espacios de alta dimensión.
  • Tiene un buen rendimiento incluso cuando algunas variables explicativas no aportan información a la predicción.
  • No produce estimaciones sesgadas ni conduce al sobreajuste.

Random Forest

\[ \hat{\theta}^B(x) = \frac{1}{B} \sum_{b=1}^B{t^*_b(x)} \] \(\hat{\theta}^B(x)\) es la predicción final

\(b\) es la muestra individual de arranque

\(B\) es el número total de árboles

\(t^*_b\) es el árbol de decisión individual

Random Forest

Support Vector Machines

Support Vector Machines

Support Vector Machines (SVM; Hearst et al., 1998, Steinwart et al., 2008) tiene como objetivo encontrar un hiperplano en un espacio de N dimensiones que clasifique distintamente los puntos de datos dados. Como hay muchos hiperplanos que se pueden seleccionar, SVM maximiza la distancia de margen entre puntos de datos de diferentes clases.

Redes Neuronales

Redes Neuronales

Las Redes Neuronales (NN), también conocidas como Redes Neuronales Artificiales, están compuestas por capas de nodos. Tienen una capa de entrada, una o más capas ocultas y una capa de salida. Cada nodo se conecta a otro y tiene un peso y un umbral asociados (IBM). Si la salida de un nodo diseñado está por encima del umbral especificado, el nodo envía información al siguiente.

Redes Neuronales

Uso de machine learning para la predicción

Uso de machine learning para la predicción

Ahora que hemos obtenido una visión general de diferentes algoritmos de machine learning, vamos a usarlos en un ejercicio de regresión. El conjunto de datos mtcars (incluido en R) comprende el consumo de combustible y 10 aspectos del diseño y rendimiento de automóviles para 32 automóviles.

print(mtcars)
##                      mpg cyl  disp  hp drat    wt  qsec vs am gear carb
## Mazda RX4           21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4
## Mazda RX4 Wag       21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4
## Datsun 710          22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1
## Hornet 4 Drive      21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1
## Hornet Sportabout   18.7   8 360.0 175 3.15 3.440 17.02  0  0    3    2
## Valiant             18.1   6 225.0 105 2.76 3.460 20.22  1  0    3    1
## Duster 360          14.3   8 360.0 245 3.21 3.570 15.84  0  0    3    4
## Merc 240D           24.4   4 146.7  62 3.69 3.190 20.00  1  0    4    2
## Merc 230            22.8   4 140.8  95 3.92 3.150 22.90  1  0    4    2
## Merc 280            19.2   6 167.6 123 3.92 3.440 18.30  1  0    4    4
## Merc 280C           17.8   6 167.6 123 3.92 3.440 18.90  1  0    4    4
## Merc 450SE          16.4   8 275.8 180 3.07 4.070 17.40  0  0    3    3
## Merc 450SL          17.3   8 275.8 180 3.07 3.730 17.60  0  0    3    3
## Merc 450SLC         15.2   8 275.8 180 3.07 3.780 18.00  0  0    3    3
## Cadillac Fleetwood  10.4   8 472.0 205 2.93 5.250 17.98  0  0    3    4
## Lincoln Continental 10.4   8 460.0 215 3.00 5.424 17.82  0  0    3    4
## Chrysler Imperial   14.7   8 440.0 230 3.23 5.345 17.42  0  0    3    4
## Fiat 128            32.4   4  78.7  66 4.08 2.200 19.47  1  1    4    1
## Honda Civic         30.4   4  75.7  52 4.93 1.615 18.52  1  1    4    2
## Toyota Corolla      33.9   4  71.1  65 4.22 1.835 19.90  1  1    4    1
## Toyota Corona       21.5   4 120.1  97 3.70 2.465 20.01  1  0    3    1
## Dodge Challenger    15.5   8 318.0 150 2.76 3.520 16.87  0  0    3    2
## AMC Javelin         15.2   8 304.0 150 3.15 3.435 17.30  0  0    3    2
## Camaro Z28          13.3   8 350.0 245 3.73 3.840 15.41  0  0    3    4
## Pontiac Firebird    19.2   8 400.0 175 3.08 3.845 17.05  0  0    3    2
## Fiat X1-9           27.3   4  79.0  66 4.08 1.935 18.90  1  1    4    1
## Porsche 914-2       26.0   4 120.3  91 4.43 2.140 16.70  0  1    5    2
## Lotus Europa        30.4   4  95.1 113 3.77 1.513 16.90  1  1    5    2
## Ford Pantera L      15.8   8 351.0 264 4.22 3.170 14.50  0  1    5    4
## Ferrari Dino        19.7   6 145.0 175 3.62 2.770 15.50  0  1    5    6
## Maserati Bora       15.0   8 301.0 335 3.54 3.570 14.60  0  1    5    8
## Volvo 142E          21.4   4 121.0 109 4.11 2.780 18.60  1  1    4    2

Uso de machine learning para la predicción

Queremos predecir la potencia (hp) de algunos autos utilizando los otros parámetros de los autos:

  • mpg: Millas por galón (US)
  • cyl: Número de cilindros
  • disp: Desplazamiento (cu.in.)
  • drat: Relación del eje trasero
  • wt: Peso (1000 lbs)
  • qsec: Tiempo de 1/4 de milla
  • vs: Motor (0 = motor en V, 1 = en linea)
  • am: Transmisión (0 = automática, 1 = manual)
  • gear: Número de marchas hacia adelante

Uso de machine learning para la predicción

En este ejercicio utilizaremos Random Forest, Support Vector Machines y Neural Networks. Por lo tanto, carguemos los paquetes correspondientes:

library(randomForest)
library(e1071)
library(neuralnet)

Uso de machine learning para la predicción

Para predecir la potencia (hp) de los autos, vamos a separar aleatoriamente los datos en una muestra de entrenamiento (para entrenar los algoritmos de machine learning) y una muestra de validación (para evaluar su rendimiento). Para garantizar la reproducibilidad, estableceremos una semilla en el código. Establecer una semilla en R significa inicializar un generador de números pseudoaleatorios.

set.seed(28)

Uso de machine learning para predicción

Para separar los datos, vamos a almacenar en un objeto el número total de autos.

n <- nrow(mtcars)

Después de esto, podemos generar una muestra aleatoria de posiciones. Vamos a utilizar el 80% del número total de autos con fines de entrenamiento.

pos_training <- sample(1:nrow(mtcars), round(n*0.8))
print(pos_training)
##  [1] 17 32  9 31  2 14 11  1 19  8  3 25 16 20 22 12 30 18  7 24 27  6 29 23 28
## [26]  5

Uso de machine learning para predicción

Ahora podemos separar el conjunto de datos mtcars en conjuntos de entrenamiento y validación.

training   <- mtcars[pos_training,]
validation <- mtcars[-pos_training,]
print(training)
##                      mpg cyl  disp  hp drat    wt  qsec vs am gear carb
## Chrysler Imperial   14.7   8 440.0 230 3.23 5.345 17.42  0  0    3    4
## Volvo 142E          21.4   4 121.0 109 4.11 2.780 18.60  1  1    4    2
## Merc 230            22.8   4 140.8  95 3.92 3.150 22.90  1  0    4    2
## Maserati Bora       15.0   8 301.0 335 3.54 3.570 14.60  0  1    5    8
## Mazda RX4 Wag       21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4
## Merc 450SLC         15.2   8 275.8 180 3.07 3.780 18.00  0  0    3    3
## Merc 280C           17.8   6 167.6 123 3.92 3.440 18.90  1  0    4    4
## Mazda RX4           21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4
## Honda Civic         30.4   4  75.7  52 4.93 1.615 18.52  1  1    4    2
## Merc 240D           24.4   4 146.7  62 3.69 3.190 20.00  1  0    4    2
## Datsun 710          22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1
## Pontiac Firebird    19.2   8 400.0 175 3.08 3.845 17.05  0  0    3    2
## Lincoln Continental 10.4   8 460.0 215 3.00 5.424 17.82  0  0    3    4
## Toyota Corolla      33.9   4  71.1  65 4.22 1.835 19.90  1  1    4    1
## Dodge Challenger    15.5   8 318.0 150 2.76 3.520 16.87  0  0    3    2
## Merc 450SE          16.4   8 275.8 180 3.07 4.070 17.40  0  0    3    3
## Ferrari Dino        19.7   6 145.0 175 3.62 2.770 15.50  0  1    5    6
## Fiat 128            32.4   4  78.7  66 4.08 2.200 19.47  1  1    4    1
## Duster 360          14.3   8 360.0 245 3.21 3.570 15.84  0  0    3    4
## Camaro Z28          13.3   8 350.0 245 3.73 3.840 15.41  0  0    3    4
## Porsche 914-2       26.0   4 120.3  91 4.43 2.140 16.70  0  1    5    2
## Valiant             18.1   6 225.0 105 2.76 3.460 20.22  1  0    3    1
## Ford Pantera L      15.8   8 351.0 264 4.22 3.170 14.50  0  1    5    4
## AMC Javelin         15.2   8 304.0 150 3.15 3.435 17.30  0  0    3    2
## Lotus Europa        30.4   4  95.1 113 3.77 1.513 16.90  1  1    5    2
## Hornet Sportabout   18.7   8 360.0 175 3.15 3.440 17.02  0  0    3    2
print(validation)
##                     mpg cyl  disp  hp drat    wt  qsec vs am gear carb
## Hornet 4 Drive     21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1
## Merc 280           19.2   6 167.6 123 3.92 3.440 18.30  1  0    4    4
## Merc 450SL         17.3   8 275.8 180 3.07 3.730 17.60  0  0    3    3
## Cadillac Fleetwood 10.4   8 472.0 205 2.93 5.250 17.98  0  0    3    4
## Toyota Corona      21.5   4 120.1  97 3.70 2.465 20.01  1  0    3    1
## Fiat X1-9          27.3   4  79.0  66 4.08 1.935 18.90  1  1    4    1

Random Forest

El paquete randomForest implementa el algoritmo de random forest de Breiman para clasificación y regresión.

randomForest(x, data, ntree = 500, mtry, nodesize, maxnodes, na.action)

x: data frame o una matriz de predictores, o una fórmula que describe el modelo a ajustar

data: data frame que contiene las variables del modelo

ntree: número de árboles a construir

mtry: número de variables muestreadas aleatoriamente como candidatos en cada división

Random Forest

randomForest(x, data, ntree = 500, mtry, nodesize, maxnodes, na.action)

nodesize: Tamaño mínimo de los nodos terminales. Si se establece un número mayor, se generarán árboles más pequeños (y, por lo tanto, se tardará menos tiempo)

maxnodes: Número máximo de nodos terminales que pueden tener los árboles del bosque. Si no se especifica, los árboles se crecerán hasta el máximo posible (sujeto a los límites de nodesize)

na.action: Una función para especificar la acción a tomar si se encuentran valores faltantes (NA)

Existen más parámetros relacionados con esta función. Por favor, escribe ?randomForest para verlos.

Random Forest

En este ejemplo, utilizaremos los valores predeterminados de la función. La fórmula que se aplicará es \(hp\sim.\). Esto significa que \(hp\) es la variable dependiente y el \(.\) significa que se utilizarán todas las variables del conjunto de datos de entrenamiento.

rf_model <- randomForest(hp ~ ., data = training)

En caso de que queramos utilizar solo parámetros específicos, podemos configurarlos de la siguiente manera:

randomForest(hp ~ mpg + cyl + disp, data = training)

Random Forest

Podemos imprimir el objeto rf_model.

print(rf_model)
## 
## Call:
##  randomForest(formula = hp ~ ., data = training) 
##                Type of random forest: regression
##                      Number of trees: 500
## No. of variables tried at each split: 3
## 
##           Mean of squared residuals: 1322.719
##                     % Var explained: 73.54

Para predecir los valores de hp, utilizaremos la función predict, que requiere el modelo y las covariables para predecir la potencia de los automóviles. Para dejarlo explícito, excluiremos la columna de hp del data.frame.

rf_prediction <- predict(rf_model, validation[,-4])
print(rf_prediction)
##     Hornet 4 Drive           Merc 280         Merc 450SL Cadillac Fleetwood 
##          123.43118          122.20713          173.63778          219.34986 
##      Toyota Corona          Fiat X1-9 
##           96.53314           74.71494

Random Forest

Podemos crear un data.frame con los valores reales (Obs) y las predicciones del Random Forest (RF).

data.frame(Obs = validation[,4], RF = rf_prediction)
##                    Obs        RF
## Hornet 4 Drive     110 123.43118
## Merc 280           123 122.20713
## Merc 450SL         180 173.63778
## Cadillac Fleetwood 205 219.34986
## Toyota Corona       97  96.53314
## Fiat X1-9           66  74.71494

Random Forest

La función varImpPlot generará un gráfico de puntos de la importancia de las variables según lo medido por un Random Forest.

varImpPlot(rf_model)

Support Vector Machines

El paquete e1071 se utiliza para entrenar una máquina de vectores de soporte para clasificación y regresión.

svm(formula, data, kernel, degree, gamma, coef0, na.action)

formula: un data.frame o una matriz de predictores, o una fórmula que describe el modelo a ajustar

data: un data.frame que contiene las variables del modelo

kernel: el núcleo utilizado en el entrenamiento y la predicción: lineal, polinómico, radial y sigmoidal

Support Vector Machines

svm(formula, data, kernel, degree, gamma, coef0, na.action)

degree: parámetro relacionado con el núcleo polynomial

gamma: parámetro necesario para todos los núcleos excepto el lineal

coef0: parámetro necesario para los núcleos de tipo polinómico y sigmoidal

Support Vector Machines

svm(formula, data, kernel, degree, gamma, coef0, na.action)

na.action: una función para especificar la acción a tomar si se encuentran valores faltantes (NA)

Existen más parámetros relacionados con esta función. Por favor, escribe ?svm para verlos.

Support Vector Machines

Al igual que RF, la fórmula que se aplicará es \(hp\sim.\). Utilizaremos diferentes núcleos.

svm_model_lin    <- svm(hp ~ ., data = training, kernel = "linear")
svm_model_poly   <- svm(hp ~ ., data = training, kernel = "polynomial")
svm_model_rad    <- svm(hp ~ ., data = training, kernel = "radial")
svm_model_sig    <- svm(hp ~ ., data = training, kernel = "sigmoid")

Ahora podemos predecir los valores de hp:

svm_prediction_lin  <- predict(svm_model_lin, validation[,-4])
svm_prediction_poly <- predict(svm_model_poly, validation[,-4])
svm_prediction_rad  <- predict(svm_model_rad, validation[,-4])
svm_prediction_sig  <- predict(svm_model_sig, validation[,-4])

Support Vector Machines

Finalmente, podemos resumir los resultados.

data.frame(Obs = validation[,4], SVM_lin = svm_prediction_lin,
           SVM_poly = svm_prediction_poly, SVM_rad = svm_prediction_rad, 
           SVM_sig = svm_prediction_sig)
##                    Obs   SVM_lin  SVM_poly   SVM_rad   SVM_sig
## Hornet 4 Drive     110 121.51574 123.47778 106.47474 120.99467
## Merc 280           123 156.98698 146.04602 125.07476 142.84838
## Merc 450SL         180 156.65475 168.71280 171.96696 160.63893
## Cadillac Fleetwood 205 236.92939 223.57702 211.27667 226.52371
## Toyota Corona       97  45.72309  80.75757  91.49229  82.29314
## Fiat X1-9           66  69.46686  77.96054  77.80608  68.92851

Neural Networks

El paquete neuralnet permite entrenar redes neuronales utilizando retropropagación, retropropagación resiliente (RPROP) con o sin retroceso de peso.

neuralnet(formula, data, hidden, threshold, stepmax, algorithm)

formula: un data.frame o una matriz de predictores, o una fórmula que describe el modelo a ajustar

data: un data.frame que contiene las variables del modelo

hidden: un vector de enteros que especifica el número de neuronas ocultas en cada capa

Neural Networks

neuralnet(formula, data, hidden, threshold, stepmax, algorithm)

threshold: un valor numérico que especifica el umbral para las derivadas parciales de la función de error como criterio de detención

stepmax: el número máximo de pasos para el entrenamiento de la red neuronal

Neural Networks

neuralnet(formula, data, hidden, threshold, stepmax, algorithm)

algorithm: una cadena de texto que contiene el tipo de algoritmo para calcular la red neuronal. Los tipos posibles son: ‘backprop’, ‘rprop+’, ‘rprop-’, ‘sag’ o ‘slr’

Existen más parámetros relacionados con esta función. Por favor, escribe ?neuralnet para verlos.

Neural Networks

La fórmula que se aplicará es \(hp\sim.\).

nn_model  <- neuralnet(hp~., data=training)

Ahora podemos predecir los valores de hp:

nn_prediction  <- predict(nn_model, validation[,-4])

Neural Networks

Finalmente, podemos resumir los resultados.

data.frame(Obs=validation[,4], NN=nn_prediction)
##                    Obs       NN
## Hornet 4 Drive     110 150.4998
## Merc 280           123 150.4998
## Merc 450SL         180 150.4998
## Cadillac Fleetwood 205 150.4998
## Toyota Corona       97 150.4998
## Fiat X1-9           66 150.4998
  • ¿Qué ocurrió?

Neural Networks

La primera razón a considerar cuando se obtienen resultados extraños con redes neuronales es la normalización. Los datos deben estar normalizados, de lo contrario, sí, el entrenamiento dará como resultado una NN sesgada que producirá el mismo resultado todo el tiempo.

Neural Networks

Normalicemos los datos.

norm_df <-rbind(training, validation)

#define Min-Max normalization function
minmax_norm <- function(x) {
  (x - min(x)) / (max(x) - min(x))
}

norm <- apply(norm_df, 2, minmax_norm)
head(norm)
##                         mpg cyl      disp        hp      drat        wt
## Chrysler Imperial 0.1829787 1.0 0.9201796 0.6289753 0.2165899 0.9798006
## Volvo 142E        0.4680851 0.0 0.1244699 0.2014134 0.6221198 0.3239581
## Merc 230          0.5276596 0.0 0.1738588 0.1519435 0.5345622 0.4185630
## Maserati Bora     0.1957447 1.0 0.5734597 1.0000000 0.3594470 0.5259524
## Mazda RX4 Wag     0.4510638 0.5 0.2217511 0.2049470 0.5253456 0.3482485
## Merc 450SLC       0.2042553 1.0 0.5106011 0.4522968 0.1428571 0.5796471
##                         qsec vs am gear      carb
## Chrysler Imperial 0.34761905  0  0  0.0 0.4285714
## Volvo 142E        0.48809524  1  1  0.5 0.1428571
## Merc 230          1.00000000  1  0  0.5 0.1428571
## Maserati Bora     0.01190476  0  1  1.0 1.0000000
## Mazda RX4 Wag     0.30000000  0  1  0.5 0.4285714
## Merc 450SLC       0.41666667  0  0  0.0 0.2857143

Neural Networks

Ahora, podemos seleccionar los datos correspondientes.

norm_training   <- norm[1:nrow(training),]
norm_validation <- norm[-c(1:nrow(training)),]

Finalmente, podemos aplicar el modelo de NN.

nn_model_log  <- neuralnet(hp~., data = norm_training, act.fct="logistic")
nn_model_tan  <- neuralnet(hp~., data = norm_training, act.fct="tanh")

Neural Networks

Podemos graficar la NN.

plot(nn_model_log)

Podemos graficar la NN de una manera más agradable.

plot(nn_model_log, col.hidden = 'darkgreen', col.hidden.synapse = 'darkgreen',
     show.weights = FALSE, information = FALSE, fill = 'lightblue')

Neural Networks

Podemos usar los modelos para predecir la variable.

nn_prediction_log  <- predict(nn_model_log, norm_validation[,-4])
nn_prediction_tan  <- predict(nn_model_tan, norm_validation[,-4])

Para visualizar los resultados:

data.frame(Obs = norm_validation[,4], NN_log = nn_prediction_log, NN_tan = nn_prediction_tan)
##                           Obs     NN_log       NN_tan
## Hornet 4 Drive     0.20494700 0.15848709  0.280470546
## Merc 280           0.25088339 0.30682250  0.338217726
## Merc 450SL         0.45229682 0.35768535  0.359986470
## Cadillac Fleetwood 0.54063604 0.61909455  0.705239030
## Toyota Corona      0.15901060 0.10892604 -0.003904399
## Fiat X1-9          0.04946996 0.08748578  0.041047591

Neural Networks

Finalmente, podemos volver a transformar los datos de hp:

hp_val_log <- nn_prediction_log * (max(norm_df$hp) - min(norm_df$hp)) + min(norm_df$hp)
hp_val_tan <- nn_prediction_tan * (max(norm_df$hp) - min(norm_df$hp)) + min(norm_df$hp)

data.frame(Obs = validation[,4], NN_log = hp_val_log, NN_tan = hp_val_tan)
##                    Obs    NN_log    NN_tan
## Hornet 4 Drive     110  96.85185 131.37316
## Merc 280           123 138.83077 147.71562
## Merc 450SL         180 153.22495 153.87617
## Cadillac Fleetwood 205 227.20376 251.58265
## Toyota Corona       97  82.82607  50.89506
## Fiat X1-9           66  76.75848  63.61647

Using machine learning for prediction

Finalmente, creemos un resumen de los resultados.

result <- data.frame(Obs = validation[,4], RF = rf_prediction, SVM_lin = svm_prediction_lin,
              SVM_poly = svm_prediction_poly, SVM_sig = svm_prediction_sig,
              NN_log = hp_val_log, NN_tan = hp_val_tan)
print(result)
##                    Obs        RF   SVM_lin  SVM_poly   SVM_sig    NN_log
## Hornet 4 Drive     110 123.43118 121.51574 123.47778 120.99467  96.85185
## Merc 280           123 122.20713 156.98698 146.04602 142.84838 138.83077
## Merc 450SL         180 173.63778 156.65475 168.71280 160.63893 153.22495
## Cadillac Fleetwood 205 219.34986 236.92939 223.57702 226.52371 227.20376
## Toyota Corona       97  96.53314  45.72309  80.75757  82.29314  82.82607
## Fiat X1-9           66  74.71494  69.46686  77.96054  68.92851  76.75848
##                       NN_tan
## Hornet 4 Drive     131.37316
## Merc 280           147.71562
## Merc 450SL         153.87617
## Cadillac Fleetwood 251.58265
## Toyota Corona       50.89506
## Fiat X1-9           63.61647

Exercise - machine learning for streamflow prediction

Machine learning for streamflow prediction

Ahora hagamos un ejemplo más práctico. Imagina que tenemos dos cuencas anidadas y tenemos algunos datos faltantes en la cuenca inferior que necesitamos. Usemos el aprendizaje automático para predecir esos valores faltantes.

Predicción de caudal

Para trabajar con la serie temporal, utilizaremos el paquete hydroTSM. Primero, vamos a leer los datos, que se almacenan como un objeto zoo.

library(hydroTSM)
q_path <- "C:/User/Data/L7_Machine_Learning_I/Streamflow_mm.zoo"
q <- read.zoo(q_path, header = TRUE)

Predicción de caudal

plot(q, main = "Caudal [mm]")

Predicción de caudal

Separemos los datos de la cuenca y la subcuenca en consecuencia.

subcatch <- q[,1]
catch    <- q[,2] 

Dado que en este caso tenemos los datos completos, crearemos un objeto adicional con el caudal de la cuenca y eliminaremos algunos años de datos. De esta manera, podremos comparar los valores medidos con las predicciones.

catch_na <- catch

Predicción de caudal

Podemos almacenar las fechas del objeto zoo en un objeto.

dates <- index(q)
head(dates)
## [1] "1990-01-01" "1990-01-02" "1990-01-03" "1990-01-04" "1990-01-05"
## [6] "1990-01-06"
tail(dates)
## [1] "2018-12-26" "2018-12-27" "2018-12-28" "2018-12-29" "2018-12-30"
## [6] "2018-12-31"

Predicción de caudal

Ahora, estableceremos los primeros seis años de catch_na como NAs. Para este propósito, crearemos un objeto pos que almacenará la posición de los datos relacionados con ‘1995-12-31’.

pos <- grep("1995-12-31", dates)
catch_na[1:pos] <- NA

Predicción de caudal

En este sentido, el objeto catch tendrá los datos completos.

plot(catch, main = "Caudal [mm]")

Predicción de caudal

Mientras que el objeto catch_na tendrá seis años de valores faltantes.

plot(catch_na, main = "Caudal [mm]")

Predicción de caudal

Particionaremos los datos que se utilizarán para fines de entrenamiento y verificación utilizando el objeto pos que creamos anteriormente.

# El conjunto de entrenamiento debe comenzar un día después de "1995-12-31"
inicio <- pos + 1

subcatch_training <- subcatch[, inicio:length(subcatch)]
catch_training    <- catch[, inicio:length(catch)]

head(subcatch_training)
## 1990-01-01 1990-01-02 1990-01-03 1990-01-04 1990-01-05 1990-01-06 
##   4.540620   4.428097   4.269241   4.203052   4.090529   4.004482
head(catch_training)
## 1990-01-01 1990-01-02 1990-01-03 1990-01-04 1990-01-05 1990-01-06 
##   2.521102   2.452117   2.392539   2.389403   2.285925   2.220075

Predicción de caudal

training <- data.frame(subcatch = subcatch_training, catch = catch_training)
head(training)
##            subcatch    catch
## 1990-01-01 4.540620 2.521102
## 1990-01-02 4.428097 2.452117
## 1990-01-03 4.269241 2.392539
## 1990-01-04 4.203052 2.389403
## 1990-01-05 4.090529 2.285925
## 1990-01-06 4.004482 2.220075
tail(training)
##            subcatch    catch
## 2018-12-26 3.263157 1.859470
## 2018-12-27 3.203586 1.803027
## 2018-12-28 3.130777 1.774806
## 2018-12-29 3.177110 1.771670
## 2018-12-30 3.276395 1.909641
## 2018-12-31 3.005016 1.762263
summary(training)
##     subcatch           catch        
##  Min.   : 0.9664   Min.   : 0.4829  
##  1st Qu.: 2.7403   1st Qu.: 1.3672  
##  Median : 4.5208   Median : 2.8190  
##  Mean   : 5.7265   Mean   : 4.1321  
##  3rd Qu.: 7.2147   3rd Qu.: 5.7070  
##  Max.   :45.2076   Max.   :47.0982  
##  NA's   :123       NA's   :7

Predicción de caudal

Durante el entrenamiento, podemos enfrentar algunos problemas si tenemos datos faltantes. Excluyamos los datos faltantes del entrenamiento.

pos_nas <- which(is.na(training$subcatch))
pos_nas <- c(pos_nas, which(is.na(training$catch)))
print(pos_nas)
##   [1] 7396 7397 7398 7399 7400 7401 7402 7403 7404 7405 7406 7407 7408 7409 7410
##  [16] 7411 7412 7413 7414 7415 7416 7417 7418 7419 7420 7421 7422 7423 7424 7425
##  [31] 7426 7427 7428 7429 7430 7431 7432 7433 7434 7435 7436 7437 7438 7439 7440
##  [46] 7441 7442 7443 7444 7445 7446 7447 7448 7449 7450 7451 7452 7453 7454 7455
##  [61] 7456 8188 8189 9442 9443 9444 9445 9446 9447 9448 9449 9450 9451 9453 9454
##  [76] 9455 9456 9457 9458 9460 9461 9462 9463 9464 9465 9466 9467 9468 9469 9470
##  [91] 9471 9472 9473 9474 9475 9476 9477 9478 9479 9480 9481 9482 9483 9484 9485
## [106] 9486 9487 9488 9489 9490 9491 9492 9493 9494 9495 9496 9497 9498 9499 9798
## [121] 9799 9800 9801 2237 2238 2239 2240 5537 6686 9831

training <- training[-pos_nas,]
summary(training)
##     subcatch           catch        
##  Min.   : 0.9664   Min.   : 0.4829  
##  1st Qu.: 2.7403   1st Qu.: 1.3609  
##  Median : 4.5274   Median : 2.8692  
##  Mean   : 5.7287   Mean   : 4.1585  
##  3rd Qu.: 7.2147   3rd Qu.: 5.7383  
##  Max.   :45.2076   Max.   :47.0982

Predicción de caudal

El siguiente paso es entrenar los modelos. Para este propósito, utilizaremos Random Forest y Support Vector Machines con los diferentes kernels.

model_rf         <- randomForest(catch ~ ., data = training)
model_svm_lin    <- svm(catch ~ ., data = training, kernel = "linear")
model_svm_poly   <- svm(catch ~ ., data = training, kernel = "polynomial")
model_svm_rad    <- svm(catch ~ ., data = training, kernel = "radial")
model_svm_sig    <- svm(catch ~ ., data = training, kernel = "sigmoid")

Predicción de caudal

Ahora, preparemos los datos que se utilizarán para predecir los valores para la subcuenca (es decir, el caudal diario de los primeros seis años de la cuenca).

covariate <- data.frame(subcatch = subcatch[1:pos])
head(covariate)
##            subcatch
## 1990-01-01 4.540620
## 1990-01-02 4.428097
## 1990-01-03 4.269241
## 1990-01-04 4.203052
## 1990-01-05 4.090529
## 1990-01-06 4.004482
tail(covariate)
##            subcatch
## 1995-12-26 3.918436
## 1995-12-27 3.819151
## 1995-12-28 3.653676
## 1995-12-29 3.541154
## 1995-12-30 3.441869
## 1995-12-31 3.375679

Predicción de caudal

Finalmente, podemos predecir los valores de caudal utilizando los diferentes modelos y el objeto covariable.

prediction_rf       <- predict(model_rf, covariate)
prediction_svm_lin  <- predict(model_svm_lin, covariate)
prediction_svm_poly <- predict(model_svm_poly, covariate)
prediction_svm_rad  <- predict(model_svm_rad, covariate)
prediction_svm_sig  <- predict(model_svm_sig, covariate)

Predicción de caudal

Ahora que tenemos los valores predichos, podemos convertirlos en objetos zoo.

prediction_rf       <- zoo(prediction_rf, dates[1:pos])
prediction_svm_lin  <- zoo(prediction_svm_lin, dates[1:pos])
prediction_svm_poly <- zoo(prediction_svm_poly, dates[1:pos])
prediction_svm_rad  <- zoo(prediction_svm_rad, dates[1:pos])
prediction_svm_sig  <- zoo(prediction_svm_sig, dates[1:pos])

Predicción de caudal

Podemos graficar los resultados de las predicciones de la siguiente manera:

plot(catch_na)
lines(prediction_rf, col="green")
lines(prediction_svm_lin, col="blue")
lines(prediction_svm_poly, col="red")
lines(prediction_svm_rad, col="cyan")
lines(prediction_svm_sig, col="purple")

grid()
legend("topright", lwd = 1, c("RF", "SVM-Lineal", "SVM-Polinomial", "SVM-Radial", "SVM-Sigmoid"),
       col = c("black", "green", "blue", "red", "cyan", "purple"), bty = "n")

Predicción de caudal

Podemos graficar los resultados de las predicciones de la siguiente manera:

Predicción de caudal

Eliminando los modelos de bajo rendimiento:

plot(catch_na)
lines(prediction_rf, col="green")
lines(prediction_svm_lin, col="blue")
lines(prediction_svm_rad, col="cyan")

grid()
legend("topright", lwd = 1, c("RF", "SVM-Lineal", "SVM-Radial"),
       col = c("black", "green", "blue", "cyan"), bty = "n")

Predicción de caudal

Podemos graficar los resultados de las predicciones de la siguiente manera:

Predicción de caudal

Podemos utilizar un indicador de rendimiento para evaluar qué modelo tuvo el mejor desempeño. La eficiencia de Kling-Gupta (KGE; Gupta et al., 2009, Kling et al., 2012) tiene tres componentes: el coeficiente de correlación de Pearson (\(r\)); la relación de sesgo (\(\beta\)); y la relación de variabilidad (\(\gamma\)).

Predicción de caudal

\[\scriptsize \textrm{KGE'} = 1 - \sqrt{ (r-1)^2 + (\beta-1)^2 + (\gamma-1)^2 } \]

 

\[\scriptsize r = \frac{ \sum_{i=1}^n{(O_i - \bar{O}) (S_i - \bar{S}) } }{\sqrt{\sum_{i=1}^n{(O_i - \bar{O})^2}}\sqrt{\sum_{i=1}^n{(S_i - \bar{S})^2}}}; \beta = \frac{\mu_{s}}{\mu_{o}}; \gamma = \frac{CV_{s}}{CV_{o}} = \frac{\sigma_{s}/\mu_{s}}{\sigma_{o}/\mu_{o}} \]

Donde \(\mu\) es el promedio del caudal (\(Q\)), \(CV\) es el coeficiente de variación, \(\sigma\) representa la desviación estándar de \(Q\), y los subíndices \(s\) y \(o\) representan el caudal simulado y observado, respectivamente. El KGE’ y sus componentes tienen su valor óptimo en uno. \(r\) se relaciona a la dinámica temporal, mientras \(\beta\) y \(\gamma\) se relacionan a el volumen y la variabilidad, respectivamente.

Predicción de caudal

Podemos crear un objeto que almacene las observaciones del período de verificación.

obs <- catch[1:pos]

Finalmente, podemos aplicar el KGE’ utilizando la función KGE del paquete hydroGOF.

library(hydroGOF)

Predicción de caudal

KGE(obs = obs, sim = prediction_rf, method = "2012", out.type = "full")
## $KGE.value
## [1] 0.9234464
## 
## $KGE.elements
##         r      Beta     Gamma 
## 0.9682697 1.0035594 0.9304229
KGE(obs = obs, sim = prediction_svm_lin, method = "2012", out.type = "full")
## $KGE.value
## [1] 0.9080331
## 
## $KGE.elements
##         r      Beta     Gamma 
## 0.9559886 0.9970294 0.9193025
KGE(obs = obs, sim = prediction_svm_poly, method = "2012", out.type = "full")
## $KGE.value
## [1] 0.08816593
## 
## $KGE.elements
##         r      Beta     Gamma 
## 0.6845051 0.7331579 1.8128343
KGE(obs = obs, sim = prediction_svm_rad, method = "2012", out.type = "full")
## $KGE.value
## [1] 0.8847093
## 
## $KGE.elements
##         r      Beta     Gamma 
## 0.9517110 0.9925498 0.8955749
KGE(obs = obs, sim = prediction_svm_sig, method = "2012", out.type = "full")
## $KGE.value
## [1] -90.43503
## 
## $KGE.elements
##           r        Beta       Gamma 
##  -0.7320031 -90.2282000  -4.8975169

¡Gracias por su atención!